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C ; ABSTRACT 

. As part of an effort to understand the origin of open clusters, we present a statistical 

analysis of the currently observed Pleiades. Starting with a photometric catalog of the 
cluster, we employ a maximum likelihood technique to determine the mass distribution 
of its members, including single stars and both components of binary systems. We 
' find that the overall binary fraction for unresolved pairs is 68%. Extrapolating to 

^ ■ include resolved systems, this fraction climbs to about 76%, significantly higher than 

the accepted field-star result. Both figures are sensitive to the cluster age, for which 
we have used the currently favored value of 125 Myr. The primary and secondary 
! masses within binaries are correlated, in the sense that their ratios are closer to unity 

I than under the hypothesis of random pairing. We map out the spatial variation of 

the cluster's projected and three-dimensional mass and number densities. Finally, we 
revisit the issue of mass segregation in the Pleiades. We find unambiguous evidence of 
■ segregation, and introduce a new method for quantifying it. 
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Subject headings: open clusters and associations: individual (Pleiades) — stars: mass 
function, statistics — binaries: general 



1. Introduction 

Open clusters, with their dense central concentrations of stars, are relatively easy to i dentify. 



Over a thousand syste ms are known, and the census is thought to be complete out to 2 kpc ([Brown 



2001 



Dias et al.ll2002l ). Because the clusters are no longer buried within interstellar gas and dust. 



their internal structure and dynamics is also more accessible than for younger groups. 

Despite these favorable circumstances, many basic questions remain unanswered. Most fun- 
damentally, how do open clusters form? All observed systems have undergone some degree of 
dynamical relaxation. Thus, the present-day distribution of stellar mass differs from the one just 
after disruption of the parent cloud. Recovering this initial configuration will clearly be of value in 
addressing the formation issue. But such reconstruction presupposes, and indeed requires, that we 
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gauge accurately the stellar content of present-day clusters. Even this simpler issue is a non-trivial 
one, as we show here. 

We consider one of the most intensively studied open clusters, the Pleiades. Again, our ultimate 
goal is to trace its evolution from the earliest, embedded phase to the present epoch. Focusing here 
on the latter, we ask the following questions: What are the actual masses of the member stars? 
Do they follow the field-star initial mass function? How many of the members are single, and how 
many are in binary pairs? Are the primary and secondary masses of binaries correlated? What is 
the overall density distribution in the cluster? What is the evidence for mass segregation, and how 
can this phenomenon be quantified? 



All of these questions have been addressed previously by others. iDeacon &: Hamblyi (120041 ) 
constructed a global mass function for the Pleiades. Their method was to assign masses based 
on the observed distribution of i?-magnitudes. A more accurate assessment should account for 
the photometric infiuence of binaries. Several studies directed specifically at binaries h ave probed 
selected regions for spectroscopic pairs (jRaboud &: Mermilliodlll998l : lBouvier et al.lll997l ). However, 
the overall binary fract i on has not been carefully assessed, de spite some preliminary attempts 
(jsteele fc Jameson|[l995l : bouvier et all 1 19971 : iMoraux et alJbood ). 



A fuller investigation of these and the other issues raised requires statistical methods; these 
should prove generally useful in characterizing stellar populations. Sectioii 2 des cribes our approach, 
which employs a regularized maximum likelihood technique (e.g. ICowenlll998l ). A similar method 
has been ap plied to other astronomical problems, including the rec onstruction of cloud shape s 
( Tassis 2007 ) and the investigation of binarity within globular clusters ( Romani Sz WeinberelE99ll ). 
Our study is the first to apply this versatile tool to young stellar groups. In doing so, we also relax 
many of the restrictive assumptions adopted by previous researchers. Section 3 presents our derived 
mass function for the Pleiades, along with our results for binarity. The density distribution, as well 
our quantification of mass segregation, are the topics of Section 4. Finally, Section 5 summarizes 
our findings, critically reexamines the binarity issue, and indicates the next steps in this continuing 
study. 



2. Method of Solution 

2.1. Stellar Mass Probability Function 

The basic problem is how to assign stellar masses to all the point-like sources believed to be 
cluster members. In many cases, the source is actually a spatially unresolved binary pair. More 
rarely, it is a triple or higher-order system; for simplicity, we ignore this possibility. The available 
observational data consists of photometry in several wavebands. Given the inevitable, random error 
associated with each photometric measurement, it is not possible to identify a unique primary and 
secondary mass for each source. Instead we adopt a statistical approach that finds the relative 



- 3 - 



probability for each source to contain specific combinations of primary and secondary masses. 



We introduce a stellar mass probability density, to be denoted $(m 



■pi 



This two-dimensional 

function is defined so that ^{nip, mg) Amp Am^ is the probability that a binary system exists with 
primary mass (hereafter expressed in solar units) in the interval nip to nip + Anip and secondary 
mass from nis to nis + Am^. Single stars are viewed here as binaries with = 0. We normalize 
the function over the full mass range: 



dnir. 



dnig ^{nip, nis 



1 . 



(1) 



Note that we integrate the secondary mass nig only to nip, its maximum value. Furthermore, we 
set the lower limit of nig to 0, in order to account for single stars. It is assumed that $ = for 
< nig < mmin- Here, the global minimum mass mmin is taken to be 0.08, the brown dwarf limit. 
We consider rrimax to be the highest mass still on the main sequence . The age of the Pleiades has 
been established from lithium dating as 1.25 x 10^ yr (Stauffer et al.lll998l ). This figure represents 
the main-sequence lifetime for a star of 4 Mq (jSiess et al.ll2000l ). which we adopt as ni^s.x- 



We examine separately the handful of stars that are ostensibly more massive than m^^^^, and 
hence on post-main-sequence tracks. For these 11 sources, we assigned approximate masses from 
the observed spectral types, and obtained data on their known unresolved b inary companions; this 



information was taken from the Bright Star Catalogue (iHoffleit et al.lll99ll ). These systems were 



then added by hand to our mass functions. Finally, we igi iore the brown dwarf population, thought 



to co mprise from 10 to 16% of the total system mass (jPinfield et al.l Il998l : ISchwartz &: Becklin 
20o3). 



The most direct method for evaluating ^{mp,nig) would be to guess its values over a discrete 



grid of nip- and m^- values. Each (m. 



nig) pair makes a certain contribution to the received flux in 
various wavebands. Thus, any guessed ^{nip,nig) yields a predicted distribution of fluxes, which 
will generally differ from that observed. One changes the guessed function until the observed flux 
distribution is most likely to be a statistical realization of the predicted one. 

Unfortunately, this straightforward approach is impractical. The basic difficulty is the mass- 
sensitivity of stellar luminosities. For secondary masses nig only modestly less than nip, the binary 
is indistinguishable photometrically from a single star having the primary mass. A 0.6 Mq main- 
sequence star, for exan iple, has a ET-band flux of 10.81 mag at the 133 pc distance of the Pleiades 
(jSoderblom et al.ll2005l ). Pairing this star with a 0.2 Mq secondary (which is not yet on the main 
sequence) only changes the flux to 10.68 mag. In summary, the function ^(mp,nig) evaluated in 
this way is unconstrained throughout much of the mp — nig plane. 
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2.2. Correlation within Binaries 

Since binaries of even modestly large primary to secondary mass ratio are difficult to recog- 
nize observationally, we need to infer their contribution indirectly, within t he context of a larger 



theoretical framework. The physical origin of binaries is far from settled (jZinnecker Mathieu 



200 ll ). There is a growing consensus, however, that most pairs form together, perhaps within a 



single dense core. The accumulating observations of protobinaries, i.e., tight pairs of Class or 



Class I infrared sources, have bolstered this view (e.g.. iHaisch et al.ll2004l ). 



If binaries indeed form simultaneously, but independently, within a single dense core, there is 
no reason to expect a strong correlation between the component masses. (Such a correlation would 
be expected if the secondaries formed within the primaries' circumstellar disks, for example.) A 
credible starting hypothesis, then, is that each component mass is selected randomly from its 
own probability distribution. If the formation mechanism of each star is identical, then these 
distributions are also the same. That is, we postulate that the true binary contribution to ^{nip, m^) 
is (j){mp) (j){ms), where the single-star probability density (j) is properly normalized: 

max 

(/>(m) dm = 1 . (2) 

in 

Of course, not all sources are unresolved binaries. Let h represent the fraction of sources that 
are. We suppose that this fraction is independent of stellar mass, provided the mass in question can 
represent either the primary or secondary of a pair. While this hypothesis is reasonable for low-mass 
stars, it surely fails for O- and early B-type st ars, which have an especially high multiplicity (e.g.. 



Mason et al.l 119981 : iGarcia &: Mermilliodll200lh . Such massive objects, however, are not present in 



clusters of the Pleiades age. 

Accepting the assumption of a global binary fraction, we have a tentative expression for the 
full stellar mass probability: 

$(mp,ms) = 2b4>{mp) (j){ms) + (1 - b) (j){mp) 5{'ms) . (3) 

Here, the first term represents true binaries, and the second single stars of mass nip. The factor 
of 2 multiplying the first term is necessary because of the restricted range of integration for mg in 
equation (1). That is, this integration effectively covers only half of the mp-nis plane. On the other 
hand, the normalization condition of equation (2) applies to both the primary and secondary star, 
and covers the full range of mass, from m^[^ to m-maxj for each component. 

We shall see below that the strict random pairing hypothesis, as expressed in equation (3), 
does not yield the optimal match between the predicted and observed distribution of magnitudes. 
The match can be improved, in the statistical sense outlined previously, if one allows for a limited 
degree of correlation between the primary and secondary masses within binaries. In other words, 
there is an apparent tendency for more massive primaries to be accompanied by secondaries that 
have a greater mass than would result from random selection. 
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A simple way to quantify this effect is to consider the extreme cases. If there were perfect cor- 
relation between primary and secondary masses, then the contribution to ^{mpjing) from binaries 
would be b(p{mp) S^nip — rris). With no correlation at all, ^{mp,nis) is given by equation (3). We 
accordingly define a correlation coefficient c, whose value lies between and 1. Our final expression 
for $(mp, nis) uses c to define a weighted average of the two extreme cases: 

^{mp,ms) = 2b{l - c) (t){mp) 4>{ms) + bc(p{mp) 6{mp - mg) + (1 - 6) 0(mp) . (4) 

Note that the last righthand term, representing the probability of the source being a single star, is 
unaffected by the degree of mass correlation within binaries. 



2.3. Maximum Likelihood Analysis 

2.3.1. From Masses to Magnitudes 

Reconstructing the stellar mass probability ^{mp,ms) requires that we evaluate the constants 
b and c, as well as the single-star probability (j){m). To deal with this continuous function, we 
divide the full mass range into discrete bins of width Amj. Integrating over each bin, we find j/j, 
the probability of a star's mass being in that narrow interval: 

Vi = <t>{m) dm . (5) 

We symbolize the full array of yj-values by the vector y, and similarly denote other arrays below. 
Our task, then, is to find optimal values not only for h and c, but also for all but one element of y. 
The normalization of is now expressed by the constraint 

i 

which sets the last y- value. 

For each choice of 6, c, and y, equation (4) tells us the relative probability of binaries being at 
any location in the mp — mg plane. After dividing the plane into discrete bins, each labeled by an 
index a, we define fi^ as the predicted number of systems associated with a small bin centered on 
an (mp, mg) pair. If /Utot is the total number of systems, i.e., of unresolved sources in all magnitude 
bins, then our chosen 6, c, and y^-values yield the relative fractions n/fHot- 

As an example, consider a bin a in which mp and mg have different values lying between mmin 
and mmax- Then the system is an unequal-mass binary, for which equation (4) gives 

= 2b{l-c)ypys. (7) 

Mtot 

Here, yp is the element of y corresponding to the selected mp, while ys is similarly associated with 
mg. For a bin where mp = mg, the corresponding relation is 

^ = b{l-c)yl + bcyp. (8) 
tHot 
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Note the additional term accounting for correlated binaries. Note also that a factor of 2 has been 
dropped from equation (4), since we are integrating only over that portion of the mass bin with 
nis < rup. Finally, if the system is a single star, so that = 0, we have 



{l-h)yp. (9) 



Our observational data consists of a catalog of ritot sources, each of which has an apparent 
magnitude in at least two broadband filters. (In practice, these will be the /- and i^-bands; see 
below.) As before, we divide this two-dimensional magnitude space into small bins. Our choice of 
6, c, and y leads not only to a predicted distribution in mass space, but also in magnitude space. 

Let be the predicted number of sources in each magnitude bin, now labeled by the index (5. 
Then we may write the transformation from the mass to the magnitude distribution as 



V(5 = ^TipaiJ'a , (10) 

a 



which may be recast in the abbreviated form 

V = TZfi . (11) 

Here, TZ is the response matrix, whose elements Tla/3 give the probability that a source in a mass 
bin a is observed in a magnitude bin /3. In detail, this probability utilizes a theoretical isochrone 
in the color- magnitude diagram (see Section 3.1). We must also account for random errors in the 
measured photometry. In other words, a given magnitude pair can have contributions from a range 
of mass pairs. It is for this reason that each element of involves a sum over all a- values. 

We previously showed how to obtain the relative mass distribution fi/fitot, not the actual /x 
itself. However, it is the latter that we need for equation (11). To find /Utot; we sum equation (10) 
over all /3-values, and demand that this sum be ntot; the total number of observed sources: 

(3 P a 

SO that 

ntot = pttot y^y^'^/3a {—] . (12) 

pa v/^toty 

In summary, choosing b, c, and y gives us fi/ntot through through equations (7), (8) and (9). We 
then solve equation (12) for //tot- Supplied with knowledge of /i, we finally use equation (11) to 
compute 



2.3.2. Likelihood and Regularization 



Having chosen 6, c, and y, how do we adjust these so that the predicted and ob served n iagni- 
tude distributions best match? Our technical treatment here closely follows that in ICowenI (|l998l . 
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Chapter 11), but specialized to our particular application. Let np be the number of sources actually 
observed in each two-dimensional magnitude bin. We first seek the probability that the full array 
n is a statistical realization of the predicted i^. The supposition is that each clement z^^ represents 
the average number of sources in the appropriate bin. This average would be attained after a large 
number of random realizations of the underlying distribution. If individual observed values follow 
a Poisson distribution about the mean, then the probability of observing sources is 

This probability is highest when is close to u/j. 

The likelihood function L is the total probability of observing the full set of n^-values: 

L^H P{np) . (14) 

We will find it more convenient to deal with a sum rather than a product. Thus, we use 

InL = ^ n^lnz/^ - i/^ - ln(ra^!) . (15) 

13 

The strategy is then to find, for a given n, that u which maximizes InL. For this purpose, we 
may neglect the third term in the sum, which does not depend on u. We thus maximize a slightly 
modified function: 

InL' = ^ nplni^p - vp . (16) 

13 

Since, for a given n, each P {np) peaks at vp = np, maximizing InL' is equivalent to setting 
v equal to n in equation (11), and then inverting the response matrix to obtain fi. Such a direct 
inversion procedure typically yields a very noisy /x, including unphysical (negative) elements. The 
solution is to regularize our result by employing an entropy term S: 

S= -Y,yiln{yi) . (17) 

i 

The function S is largest when the elements of y are evenly spread out. Adding S to InL' and 
maximizing the total guarantees that the y- values are smoothly distributed, i.e., that 4>{m) is also 
a smooth function. 

In practice, we also want to vary the relative weighting of S and InL'. We do this by defining 
a regularization parameter A, and then maximizing the function F, where 

F = AlnL' + S . (18) 

For any given value of A, maximizing F yields an acceptably smooth y that reproduces well the 
observed data. For the optimal solution, we find that value of A which gives the best balance 
between smoothness of the derived 0(m) and accuracy of fit. We do this by considering another 
statistical measure, the bias. 
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2.3.3. Minimizing the Bias 

Our observational dataset, n, is an imperfect representation of the unknown probability density 
4>{m) in two senses. As already noted, n may be regarded as only one particular realization of the 
underlying distribution. Even this single realization would directly reveal cj){m) (or, equivalently, 
y) if the sample size were infinite, which of course it is not. 

Imagine that there were other realizations of (j){m). For each, we employ our maximum likeli- 
hood technique to obtain y. Averaging each yi over many trials yields its expectation value, E{yi). 
However, because of the finite sample size, E{yi) does not necessarily converge to the true value, 
ytrue^ Their difference is the bias, bi'. 

h ^ E{y,) - yf- . (19) 



The values of the biases, collecti vely denoted b, reflect the sensitivity of the estimated y to 
changes in n. Following ICowenI (|l998l . Section 11.6), we define a matrix C with elements 

The bias is then given by 

b = C{u - n) . (21) 

To evaluate the derivatives in C, we consider variations of the function F about its maximum. In 
matrix notation, 

C = -A-^B , (22) 

where the matrix A has elements 

S^F 

a^F 

dyi dnp 

Since F is a known function of both the y and n, the derivatives appearing in both A and B 
may be evaluated analytically. Another matrix that will be useful shortly is P, whose elements 

= ^ (25) 
dyi 

are also known analytically from equations (7)- (9). 

To determine the regularization parameter A appearing in F, we seek to minimize the biases. 
In practice, we consider the weighted sum of their squared values: 



and B is given by 
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and vary A to reduce this quantity^ Here, the Wa are diagonal elements of W, the covariance of 
the biases. Recall that the elements of b are here considered to be random variables that change 
with different realizations. 

We find the covariance matrix W by repeated application of the rule for error propagation 
(|Cowenlll998l . Section 1.6). We begin with V, the covariance of n. Since these values are assumed 
to be independently, Poisson-distributed variables, V has elements 

Va/3 = vpSctp . (27) 

Here we have used the fact that the variance of the Poisson distribution equals its mean, vp. In 
this equation only, both a and (3 range over the i/- values, i.e., V is a square matrix. 

We next consider 3^, the covariance of y. This is given by 

y = CVC^ . (28) 

Finally, we obtain the desired W by 

W = TyT'^ . (29) 
The matrix in this last equation has elements 

- 1^ • (30) 
Differentiating equation (21) and applying the chain rule, we find J- to be 

T = CnV-I, (31) 

where I is the identity matrix. 



2.4. Calculation of Radial Structure 

Thus far, we have focused on determining global properties of the cluster, especially the mass 
function (f)(m). We also want to investigate the spatial distribution of stellar masses. For this 
purpose, we need not perform another maximum likelihood analysis. The reason is that we can 
treat the mass distribution at each radius as a modification of the global result. 

We divide the (projected) cluster into circular annuli, each centered on a radius r. What is fi^, 
the number of sources in an annulus that are within mass bin a? (As before, each bin is labeled 
by the masses of both binary components.) The quantity we seek is 

^^a = Y.Q-P^P- (32) 
/9 



^In practice, we require that xl be reduced to A'^, the number of free parameters in our fit. As noted by'Cowery 
( 19981 . Section 11.7), the average 6i-value at this point is about equal to its standard deviation, and so is statistically 
indistinguishable from zero. 
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Here, Q^/? is the probability that a source observed within magnitude bin (3 has component masses 
within bin a. In principle, this probability depends on radius. For example, the source could have 
a high probability of having a certain mass if there exist many such stars in that annulus, even 
stars whose real magnitude is far from the observed one. If we discount such extreme variations in 
the underlying stellar population, then we may approximate Q by its global average. 

The factor i/^ in equation (32) is the estimated number of sources at radius r in magnitude bin 
/3. We only compute, via the maximum likelihood analysis, z^^, the estimated number of sources 
in the entire cluster. But we also know n^, the observed number of sources in the annulus. If the 
total number of observed sources, n^, is non-zero, then we take 

.5.^.,. (33) 

In case np = 0, then rf^ vanishes at all radii. We then assume 

.l^=^.p, (34) 
'Hot 

where n^^^ is the observed source number of all magnitudes in the annulus. In the end, equation (32) 
attributes the radial mass variation to changes in the local magnitude distribution, rather than to 
improbable observations of special objects. 

It is clear that the global Q must be closely related to the response matrix TZ, which is the 
probability that a source with a given mass has a certain magnitude. The precise relation between 
the two follows from Bayes Theorem: 

Q /iQ,//itot /or\ 
a/3 = I^Pa 1 • (35) 

The numerator of the fraction is the probability that a source at any radius lies within the mass 
bin a, while the denominator is the probability of it lying within magnitude bin j3. Note that //tot 
and Viot, need not be identical. The first quantity is the estimated number of sources covering all 
possible masses. The second is the observed number in the magnitude range under consideration. 
In practice, this range is extensive enough that the two are nearly the same. We thus write 

Qa/3 = ^/3a — . (36) 

Using this last equation, along with equation (33), equation (32) now becomes 

T 

= Ma . (37) 

For those terms where np = 0, equation (34) tells us to replace the ratio n^/np by nl^^/ritot- 
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Summing fi"^ over all a and dividing by the area of the annulus gives the projected surface 
number density of sources as a function of radius. The total mass in each annulus is 



Am'' 



E 



(38) 



where rria is the sum of the masses of both binary components in the appropriate bin. Division of 
Am'" by the annulus area gives the projected surface mass density. Under the assumption of spher- 
ical symmetry, the c orresponding volume densi ties are then found by the standard transformation 
of the Abel integral (jBinney Merrifieldl[l99a . Section 4.2.3). 



3. Application to the Pleiades: Global Results 
3.1. The Response Matrix 



We begin with the observational data. Figure 1 is a dered dened (/, I — K) colo r-magnitude 
diagram for the Pleiades, taken from the recent compilation bv IStauffer et al.l (|2007l ). Shown are 
all sources which have high membership probabi li ty, as gauged by their colors, radial velocities, 
and proper motions (see, e.g., iDeacon &: Hamblvl (|2004l ) for one such proper motion study.) The 
lower open circles correspond to probable brown dwarfs; we exclude such objects from our study. 
Most brown dwarfs are too faint to be observed, and the population, in any case, is more sparsely 
sampled. (The magnitude cutoffs corresponding to a 0.08 Mq object are Mj = 12 and Mk = 9.) 
After also exluding the 11 bright, post-main-sequence stars, shown here as large, filled circles, we 
have a total sample size of ntot = 1245. 

The solid curve near the lower boundary of the stellar dis tribution is a combination of the 
theoretical zero-age main sequence for m^, > 1 (jSiess et al.ll2000l ) and, for lower-mass stars, a pre- 
main-sequence isochrone ( Baraffe et al 19981 ) Fl The isochrone is that for the measured cluster age 
of 125 Myr. Our basic assumption is that the observed scatter about this curve stems from two 
effects - binarity and intrinsic errors in the photometric measurements. We do not consider, for 
now, possible uncertainty in the cluster's age. (See Section 5 for the effect of this uncertainty.) We 
also ignore the finite duration of star formation. This duration is roug hly 10^ yr (IPalla fc Stahleii 
2OO0l ). or about 10% of the cluster agei 



^Both theoretical results are presented in magnitudes. We have applied corrections to the the oretical _ft '-band 
magnitudes to make them consistent with the 2MASS isTs-band used in Stauffer's catalog. See iCohen et al. I (200a ) 
for this transformation. 



^We further ignore the effect of differential reddening across the cluster. IStauffer et al.l ( 20071 ) adjusted individually 
the fluxes from sources in especially obscured regions, bringing their effective extinction to the observed average Ay 
of 0.12. We therefore constructed Figure 1 by applying uniformly the corresponding Aj- and values of 0.06 and 
0.01, respectively. 
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After do i ng a polynomial fit to the mass-magnitude relations found bv lSiess et alj (|200d ) and 
BarafFe et all (j 19981 ). we have analytic expressions for Mj{mp,ms) and M^{Tnp,ms), the absolute 
/- and ii'-magnitudes for a binary consisting of a primary mass nip and secondary mg. Here, 
the superscripts indicate that the magnitudes are theoretically derived. Both Mj and are 
calculated by appropriately combining the individual absolute magnitudes for nip and m^. 



We do not actually observe Mj or for any source. What we have are dereddened, apparent 
magnitudes in these wavebands. Using the Pleiades distance of 133 pc, these apparent magnitudes 
are readily converted to absolute ones, Mj and Mk- The salient question is: Given a source with 
intrinsic magnitudes Mj and (or, equivalently, with masses nip and nig), what is the probability 
that it is observed to have magnitudes Mj and Mk'? 

Here we confront the issue of photometric errors. We assume the errors in the two wavebands 
to be normally distributed. Then the relevant probability density is 



S {Mi,MK;mp,nis 



2 TT fj/ (Tx 



exp 



{Mi - Mjy (M, 



K 



M 



Kl 



2aj 



1a\ 



(39) 



Here, S /S.Mi AMk is the probability of observing a source in magnitude bin (3, centered on 
{Mj,Mk), and having widths AM/ and AMk- 

The quantities aj and a x in equation (39) are the standard deviations of the photometric 
measurements. According to IStauffer et al.l (|2007l ). the average standa r d dev iation in the /-band 
is about 0.15. Figure 2, constructed from Table 2 of IStauffer et alJ (|2007l ). shows that ax is 
generally lower, and rises steeply with Mk for the dimmest sources The two branches of the 
curve presumably represent the results from two different observations. We do a polynomial fit to 
the upper, majority, branch, and thus have an explicit expression for axiMx)- 

Suppose now that nip and nig are centered within a mass bin a, which has widths Arup 
and Anig. Then the response matrix TZa/s is obtained by integrating S{Mj , Mx',nip,nis) over the 
magnitude bin, then averaging over the mass bin: 



n. 



■af3 



Anip Anig 



(40) 



The magnitude integrals can be expressed in terms of error functions if we reinterpret uk as being a 
function of rather than Mk- The remaining numerical integrals over nip and are performed 
by finding, for each {nip,nis) pair, the Mj- and M^-values from our polynomial fits to the mass- 
magnitude relations. 



^This rise in ax occurs because the observed is'-magnitudes are approaching the sensitivity limit of the observa - 
tions. Many of the J-band measurements come from POSS II plates , for w hich the limit is 18.5 (|Hamblv et al.ll 19931 ). 



Another large source of data was the observations of Pinfield et al. ( 2000l ) . whose limiting magnitude was 19.7. Our 



lower cutoff for brown dwarfs corresponds to an apparent /-magnitude of 17.7, so the rise in our ai should be modest. 
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3.2. Summary of Procedure and Synthetic Data Tests 

With the response matrix in hand, we are ready to input the catalog of source magnitudes. 
Before turning to the Pleiades itself, we first employed a number of synthetic datasets, in order to 
test various aspects of the code. We shall describe these tests shortly. First, however, we summarize 
the standard procedure we adopted for the analysis of any cluster, real or synthetic. 

The basic problem, we recall, is to guess the optimal values of b, c, and y that maximize the 
function T, as given in equation (18). The entropy part of F, labeled S, is directly a function of y 
(eq. (17)), while the modified likelihood function L' depends on the observed magnitude distribution 
n and the guessed one v (see eq. (16). The guessed y does not yield v itself, but the guessed mass 
distribution fi, through equations (7)-(9). It is in this transformation that the binary fraction b and 
correlation coefficient c appear. Finally, v is obtained from fi via the response matrix (eq. (11)). 

We begin the maximization procedure by first setting the regularization parameter A to zero. 
Since F = 5" in this case, the optimal set of ^/-values will be uniformly distributed, while subject to 
the normalization constraint of equation (6). We guess b, c, and y, and vary them to ma ximize F. 



For the actual maximization, we employ a standard simplex algorithm (jPress et al.l 120021 . Chapter 
10). The resulting best-fit parameters are then perturbed and the maximization rerun. This check, 
which may be redone several times, is done both to confirm convergence and to avoid becoming 
trapped in small, local maxima of the function F. We record the relevant covariances and biases, 
to be used in estimating errors in predicted quantities and to set the optimal A-value. 

The next step is to increase A slightly. We maximize F in the same way as before, again 
recording covariances and biases. We again increase A, repeating the entire procedure until xii the 
weighted sum of the biases, starts to become acceptably small. At this point, the best-fit 6, c, and 
y have been established. 

As a first test of the procedure, we introduced an artificial cluster whose single-star probability, 
(f){m), we selected beforehand. Sources were chosen randomly to have masses according to this 
distribution. In a certain fraction of the sources, our preset binary fraction b, a second star was 
added to the source. The mass of this object was also randomly chosen from (j){m). Thus, the 
correlation coefficient c was initially zero. Given both masses in a source, its intrinsic and 

are readily obtained. These magnitudes are smeared out into neighboring bins according 
to Gaussian errors with the appropriate standard deviations aj and ax- Thus, the "observed" 
magnitude distribution, n, is established. 

Figure 3 shows two representative examples. In the left panel, the chosen (^{m), is a power 
law: (p{m) oc m~^'^. On the right, we used a log-normal distribution: 

(logm — logmo)^"' 



(m) = — exp 
m 



(41) 



where C is the normalization constant. The central mass was chosen as nio = 0.2 and the loga- 
rithmic width was 0.4. The binary fraction b was chosen to be 0.30 in the power-law example, 
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and 0.68 for the log-normal distribution. The total source number ritot was 10,000 in both cases. 

The smooth curve in both panels is 0(?tt,), while the data points are the best-fit values of 
Hi/Ami, where Ami is the bin width. Shown also is the estimated error for each value. This was 
derived from the covariance matrix y, introduced in Section (2.1). Specifically, the plotted error 
is y/yii/Arrii. We divide each i/i and its associated error by Ami because 7/j is integrated over the 
bin (eq. (5)). 

It is evident that the code reproduces well the assumed (f)(m) in these two examples. Note 
that most of the scatter seen in both plots, especially in the left panel, was already present in the 
input data, which were finite realizations of the analytic distributions. The derived (i.e., predicted) 
binary fractions, b = 0.293 it 0.008 and b = 0.672 it 0.011, respectively, are also in good agreement. 
We had similar success when we reduced ntot to 1245, the actual number in our Pleiades source 
catalog. In this smaller sample, the errors in our predicted mass function and binary fraction 
increased, roughly as ntot~^^'^- 

Figure 4, taken from a dataset with ntot = 1245, shows in more detail how the regularization 
parameter A was chosen. The figure also illustrates some of the subtlety involved in this procedure. 
Plotted here, as a function of A, is xl^ defined in equation (26). As A is gradually increased, xl 
takes a sudden, sharp dip. After climbing back, xl then more slowly declines, eventually falling 
below = 21, the number of tunable parameters in this maximization (b, c, and 19 values). 

It is the second threshold (A = 0.027 in this case) that marks the true stopping point. The 
earlier dip in Xb is due, not to a decrease in the biases, but to a sharp increase in the covariances 
yV. This increase commonly occurs when the likelihood term In L' starts to become comparable to 
the entropy S in the full function T. At that point, y makes an abrupt shift away from its earlier, 
nearly uniform, distribution. With further increase in A, y settles down gradually to its optimal 
form. 

Continuing our synthetic data tests, we next introduced a correlation between the primary and 
se condary masses. First , we g enerated uncorrelated pairs, as above. Generalizing the prescription 



of Ide La Fuente Marcos! (|l996l ). we then altered the secondary mass in each source according to 

^ m J ^ . (42) 



Here, 7, a preset number between and 1, represents our imposed degree of correlation. Thus, 
setting 7 = yields the previous, uncorrelated case, while, for 7 = 1, every binary has equal-mass 
components. We ran our routine with a variety of input single-star mass functions, binary fractions, 
and degrees of correlation. 

Our general result was that the predicted y still reproduced well the synthetic (j){m). The 
binary fraction b was similarly accurate. Most significantly, the predicted c-value tracked the 
input quantity 7. Figure 5 shows this relation. We conclude that our statistical model, while 
crudely accounting for correlation by inserting a fraction of equal-mass pairs, nevertheless mimics 
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a smoother correlation, such as would be found naturally. The shaded patch in the figure is the 
probable region occupied by the real Pleiades; Secfon 3.4 below justifies this assessment!! 

One price we paid for our simplified account of correlation was that our matching of (p{m) 
was less accurate than for randomly paired input binaries. Consider, for example, the log-normal 
function of equation (41). While our best-fit y still reproduced (p{m) reasonably well, the output 
function peaked at too high a mass compared to irio- The filled circles in Figure 6 shows that 
this shift, Amo, increased with the input 7- value. Concurrently, our output function was too 
narrow compared to the input am- The (negative) difference, AcXm, displayed as open circles in 
Figure 6, was also more pronounced at higher 7. These systematic errors need to be considered 
when analyzing a real cluster. The two patched areas in the figure again represent the likely regime 
of the Pleiades, as we explain shortly. 



3.3. Empirical Mass Distributions 

We now present the results of applying our maximum likelihood a nalysis to the P l eiades itself. 



i.e., to the /- and /C-magnitudes of 1245 sources from the catalog of IStauffer et al.l (|2007l ). Our 
best-fit binary fraction was b = 0.68 it 0.02, while the correlation coefficient was c = 0.36 it 0.06. 
(These and other uncertainties represent only random statistical error, and do not include sys- 
tematic eff'ects; see Section 5.) We will discuss the implications of these findings in the following 
section. First, we examine the global distribution of stellar mass. 

The data points in Figure 7 are the best-fit values of each yi/Arrii. As in Figure 3, these points 
are a discrete representation of the single-star mass function (j){m). The large error bars on the two 
points at highest mass are due to the small number of sources gauged to be in the respective bins. 
The smooth, solid curve in Figure 7 is a log-normal mass function that best matches the empirical 
y. Referring again to equation (41), we find that nio = 0.20 it 0.04 and (Tm = 0.38 it 0.02. The 
presence of a finite binary correlation affects both estimates. Judging from Figure 6, our mo is 
overestimated by about 0.06, while am should be raised by 0.08. 

Each of our mass bins has contributions from both the primary and secondary components of 
binary pairs, as quantified by equation (4). Integrating the full stellar mass probablity $(mp,ms) 
over all secondary masses, we obtain (f)p{mp), the probability distribution of primary masses: 

4>p{mp) = / dms^{mp,ms) . (43) 
Jo 

Note that this distribution includes the possibility that the star is single (m^ = 0). 



^The prescription for mass correlation given in equation (42) is more realistic than introducing a subpopulation 
of identical-mass binaries (eq. (4)). We employed the latter device only for convenience. If we had parametrized the 
correlation through 7, equations (7) and (8) would have been numerical integrals, and the derivative matrix D in 
equation (25) would also have required numerical evaluation. 
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The solid curve in Figure 8 is a log- normal fit to the empirical (j)p{mp). Shown for comparison 
as a dashed curve is the fit for (j){m) from Figure 7. Relative to the latter function, the primary 
distribution falls off at lower masses. This falloff simply reflects the fact that less massive objects 
are more likely to be part of a binary containing a higher-mass star, and thus to be labeled as 



"secondaries. "|j In any event, we now see why the peak of (t)p{mp), rrio = 0.27 ± 0.02, is elevated 
with respect to the peak of <j)(m). Similarly, the primary distribution is also slightly narrower than 
the single-star mass function, with am = 0.35 it 0.01. 

The para r neters of our log-normal approximation to 4'p{mp) may be compared to those of 
Moraux et al.l (|2004l ). These authors fit the entire mass function. Since, however, they did not 
account for binarity, their results are more closely analogous to our primary distribution. Their 
best-fit m-o of 0.25 is close to ours, while their am of 0.52 is higher, mostly because of their iri c lusion 
of the highest-mass members. These parameters are also close to those given by IChabried ()2003l ) 
in his log-normal fit to the field-star initial mass function {mo = 0.22, a^ = 0.57). 

In Figure 9, we compare our single-star distribution 0(m) to the field-star initial mass function 
(dash ed curve). The latter, which has been raised in the figure for clarity, is taken from iKroupa 
(j200ll ). who did correct for binarity. It is apparent that (j){m) itself veers away from the IMF for 
both low- and high-mass objects. When these are added in, the resembla nce improves . The open 
circles in Figure 9 are Pleiades low-mass stars and brown dwarfs found by lBihain et all (|2006l ). We 
have normalized their data, taken from a limited area of the cluster, so that their total number of 
stars matches ours within the overlapping mass range. No s uch normalization wa s necessary for the 
11 B-type stars (filled circles), which are from the catalog of lStauffer et al.l (120071 ) but not included 
in our maximum likelihood analysis. Adding both these groups not only improves the match to the 
IMF, but also reveals a gap in the ste llar distributio r i betw een about 2 and 5 Mq. A similar gap is 
seen in the Pleiades mass function of iMoraux et al.l (|2004l . see their Figure 1). 



Our estimate for the total cluster mass, based solely on the 1245 catalog sources, is 738 Mq, 
with a 4% uncertainty. Adding in the brightest stars brings this total to 820 Mq, with the same 
relative error. Tests with synthetic data indicate that the systematic bias due to binary correlation 
raises this figure by roughly 50 Mq, to 870 Mq. Addition of the brown dwarfs would cause a 



further, relatively small, increase. For comparison, iPinfield et al.l (119981) found 735 M(^ in stars, 
and an upper limit of 131 Mq for the brown dwarf contribution. iRaboud fc MermilliodI (| 19981 ) used 
the virial theorem to estimate a total mass of 720 Mq, with a 28% uncertainty. Direct integration 
of their mass function gave 950 Mq, with an 18% fractional error. 



®We may similarly calculate a secondary mass distribution (f>s{ms) by integrating ^{mp,ms) over rrip, from rus to 
rrimax. The function 4>s has an excess of low-mass stars and drops very steeply at high masses, as most such objects 
are primaries. 



-17- 



3.4. Binarity 



The global binary fraction, b = 0.68, obtained in our analysis represents most, but not all, of 
the full binary population. Omitted here are spatially re solved systems. For th ese, the primary and 
secondary appear as separate sources in the catalog of IStauffer et al.l (|l998l ). Counting resolved 
pairs raises the total fraction to about 76%, as we now show. 

The smallest angular separation between stars in the catalog is 10". At the Pleiades distance 
of 133 pc, the corresponding physical separation is 1400 AU. An edge-on circular binary of exactly 
this orbital diameter will still be unresolved, since the components spend most of their time closer 
together. The true minimum separation in this case is 2200 AU. Here, we have divided 1400 AU 
by 2/tt, which is the average of |sin0|, for 9 randomly distributed between and 27r. 

Of course, only a relatively small fraction of binaries have separations exceeding 2200 AU. 
The average total mass of our unresolved systems is 0.71 Mq. A binary of that total mass and a 
2200 AU diameter has a period of 1.2 x 10^ yr. What fraction of binaries have eve n longer periods? 
Our a verage primary mass is 0.46 Mq, corresponding to a spectral type of Ml. iFischer &: Marcv 
(jl992l ) studied the period distribution of binaries containi ng M-type primaries. They c laimed that 
this distribution was indistinguishable from that found by iDuquennov &: Mavoii (|l99ll ) for G-type 
primaries. In this latter sample, 11% of the systems had periods greater than our limiting value. 
If the Pleiades periods are similarly distributed, then the total fraction of binaries - both resolved 
and unresolved - becomes 0.66/(1 — 0.11) = 0.76. 

Even without this augmentation, our total bina ry fraction app e ars to be inconsistent with the 
available direct observations of the Pleiades. Thus, iBouvier et al.l (119971) found visual pairs with 
periods between 40 and 3.4 x 10^ yr. Using the period distribution of iDuauennoy &: Mavoii (Il99lh 



Mermilliod et al 



to ext rapolate their observed binary fraction of 28% yields a total fraction of 60%. 
(jl992l ) observed spectroscopic pairs with periods under 3 yr. A similar exercise again yields ou7o 
We note, however, that this oste nsible concurre n ce of results is based on very broad extrapolations 
from limited data. (See Fig. 4 of lBouvier et al.l (119971 ).) 



Our derived binary fraction also exceeds that found in the field-star population. IDuquennov &: Mavoi 



(jl99ll ) found that 57% of G stars are the primaries of binary or higher-order systems. Note that 
our b represents the total probability that a star is in a binary, whethe r as the primary or sec o ndary 
component. Since G stars are rarely secondaries, the comparison with IDuquennov &: Mavoii (|l99ll ) 
is appropriate. On the other hand, M stars are frequen tly se c ondar ies, so we would expect the 
fraction of binaries with M-type primaries to be reduced. iLadal (|2006l ) has found that only 25% of 
M-stars are the primary components of binaries. Our own analysis yields a binary fraction of 45% 
for M-star primaries, still in excess of the field-star result. 

If our finding of a relatively high binary fraction proves robust, it may provide a clue to the 
progenitor state of the Pleiades and other open clusters. A similar statement applies to the corre- 
lation between component masses within binaries. Our adopted method of gauging this correlation 
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- inserting a fraction of equal-mass pairs in the mass function - is admittedly crude. Nevertheless, 
the strong result (c = 0.36 it 0.06) is significant. Referring back to Figure 5, we find that the 
Pleiades correlation is equivalent to setting 7 equal to about 0.65 in the alternative description of 
equation (42). Whatever the origin of the Pleiades binaries, the primaries and secondaries were 
not formed by completely independent processes. 

4. Application to the Pleiades: Radial Distributions 
4.1. Number and Mass Profiles 



We now employ the procedure outlined in Section 2.4 to investigate both the surface and 
volumetric density as a function of the projected distance from the cluster center. The filled circles 
in Figure 10, along with the associated error bars, represent the surface number density of so urces, 
measu red in pc"'^. The solid curve is a density profile using the empirical prescri ption of King 



2001 



1962 ). Here, the core radius is 2.1 pc, while the tidal radius is 19 pc. For comparison. lAdams et al 



also fit the surface number density profile of their low-mass stars to a King model, with a 
core radius of 2.3-3.0 pc. Our profile is also at least rou ghly consistent with the cumulative number 
distribution displayed by iRaboud h MermilliodI () 19981 ). Our best-fit tidal radius is slightly larger 



than the 17 pc cited by these authors. 

The surface mass density is plotted in an analogous fashion, again as a function of the projected 
radius. We show both the data points (small open circles) and, as the dashed curve, the best-fit 
King model. Here, the core radius is 1.3 pc, and the tidal radius is 18 pc. Note that the mass 
density profile falls off more steeply than the number density. Thus, stars near the center are 
abnormally massive, a trend we shall explore more extensively below. 

Figure 11 displays the corresponding volumetric densities. As we indicated, the deconvolution 
from surface profiles assumes sp herical symmetry. In fact, th e Pleiades is slightly asymmetric, with 
a projected axis ratio of 1.2:1 (jRaboud &: MermilliodI Il998l ). This elliptic ity is thought to stem 
from the tidal component of the Galactic gravitational field (jWielenl 1 19741 ). Under the spherical 
assumption, the filled circles and solid, smooth curve show the number density. Here, the King 
model is the same used for the surface number density in Figure 10, but deprojected into three- 
dimensional space. 

Figure 11 also shows, as the small open circles and dashed curve, the mass density as a function 
of spherical radius. Again, the King model here is the deprojected version of that from Figure 10. 
The relatively rapid falloff in the mass, as opposed to the number, density is another sign of the 
tendency for more massive stars to crowd toward the center. 

The information we used in obtaining these profiles also gives us the spatial variation of the 
binary fraction h. That is, we first used equation (37) to obtain /x*", the predicted mass distribution 
in each radial annulus. Recall that the distribution refers to both primaries and secondaries, as 
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well as single stars. The binary fraction can thus be computed locally. To within our uncertainty, 
about ±0.05 at each radial bin, we find no variation of b across the cluster. 



4.2. Mass Segregation 

We have mentioned, in a qualitative manner, that more massive cluster members tend to reside 
nearer the center. In Figure 12, we explicitly show this trend. Here, we plot {nip + nis) , the average 
system mass (primary plus secondary) as a function of the projected cluster radius. It is apparent 
that {nip + rris) monotonically falls out to about 4 pc. Beyond that point, the average mass is 
roughly constant. 

The pattern here is consistent with mass segregation, but is not a clear demonstration of that 
effect. The problem is that Figure 12 gives no indication of the relative populations at different 
annuli. If the outer ones are occupied by only a small fraction of the cluster, is mass segregation 
present? To gauge any variation in the mass distribution of stars, that distribution must be 
calculated over an adequate sample size. 



Previous authors have also claimed evidence of mass segregation, using various criteria. lAdams et al 



(j200ll ) looked a t the distribution of surface a nd volumetric number densities for a number of differ- 



ent mass bins. iRaboud MermilliodI (119981 ) divided the population by magnitude into relatively 



bright and faint stars. They calculated the cumulative number as a functi on of radius for both 



groups, and found the bright stars to be more centrally concentrated. Finally. iPinfield et al.l (|l998l ) 
fit King profiles to the surface density of various mass bins. As the average mass increases, the core 
radius shrinks. 

Figure 13 gives a simpler and more clear-cut demonstration of the effect. Here, we consider 
/at, the number of sources enclused in a given projected radius, divided by the total number of 
sources in the cluster. We also consider /m, the analogous fractional mass inside any projected 
radius. The figure then plots /m versus /at. In the absence of mass segregation, /m would equal 
/at at each annulus. This hypothetical situation is illustrated by the dotted diagonal line. 

In reality, Jm rises above /at, before they both reach unity at the cluster boundary (see upper 
smooth curve, along with the data points). This rise, as already noted, indicates that the innermost 
portion of the cluster has an anomalously large average mass, i.e., that mass segregation is present. 
Moreover, the area between the solid and dotted curves is a direct measure of the effect. In the case 
of "perfect" mass segregation, a few stars near the center would contain virtually all the cluster 
mass. The solid curve would trace the upper rectangular boundary of the plot, and the enclosed 
area would be 0.5. We thus define the Gini coefficient, G, as twice the area between the actual 
/m — fN curve and the central diagonallll For the Pleiades, we find that G = 0.20 ± 0.06. 



^The name derives from economics, where the coefficient is used to measure inequahty in the distribution of wealth 
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It is possible, at least in principle, that this effect is due entirely to a few exceptionally massive 
stars located near the center. In fact, this is not the case. We have artificially removed the 11 
brightest sources (all late-B stars) and recalculated /m versus fjy. The result is shown by the 
dashed curve in Figure 13. While the rise above the diagonal is diminished, it is still present. That 
is, the intermediate-mass population exhibits segregation, as well. 

An interesting contrast is presented by another populous group, the Orion Nebula Cluster 
(ONC). The d i stribu tion of stellar masses in this far younger system was recently studied by 
Huff &: Stahlerl (|2006l ). Figure 2 in that paper compares the stellar populations in the inner and 
outer halves of the clusterO Apart from a few high-mass objects, the two populations are essential 
identical. 

We may also construct an /m — Jn curve for the ONC, as shown here in Figure 14. The solid 
curve again lies well above the fiducial diagonal, ostensibly indicating mass segregation. However, 
removal of just the four Trapezium stars gives a dramatically different result (dashed curve) that 
is virtually indistinquishable from the diagonal. All stars except this tiny subset are similarly 
distributed. The cluster is too young to have undergone true mass segregation, a conclusion drawn 
previously from N-body simulations ([Bonnell &: Davieslll998l ). The Trapezium represents a special 
population, one that probably formed just prior to cloud dispersal (IHuff Stahleii 120071 ). 



5. Discussion 



In this paper, we have applied a versatile statistical tool, the maximum likelihood technique, to 
assess the distribution of stellar mass and the incidence of binaries in the Pleiades. We began with 
a near-infrared catalog of cluster members. Our basic assumption was that all cluster members 
share the same evolutionary age, and that any dispersion in the color-magnitude plane stems from 
binarity and random photometric errors. We were then able to infer the most probable distribution 
of masses, both for the cluster as a whole, and as a function of distance from its center. Finally, 
we introduced a simple method for gauging the degree of mass segregation in the cluster. 

One of our surprising results is the relatively high fraction of binaries. We estimate that 68% 
of all systems in the cluster are unresolved binaries; this figure climbs to about 76% if resolved 
pairs are included. These fractions are significantly higher than the accepted field-star result, so 
we should scrutinize them carefully. Could they stem from an underestimate of the photometric 
error at faint magnitudes? Since the error in / is greater than K, we artificially increased the 
dispersion aj. We kept aj at 0.15 until Mj = 9.5, below which we increased it linearly, reaching 
a I = 0.20 at Mj = 12. After redoing the maximum likelihood analysis, the global binary fraction 



( Sen 1119971 ). As we show in the Appendix, G is also half the mean mass diflFerence of radial shells, where that mass 
difference is normalized to the average system mass in the cluster. 

®Note that the axis labels in Figure 2 of Huff fc Stahlei 1 20061 ) were inadvertantly switched. 
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b for unresolved pairs is unchanged. 



Another potential difficulty is our neglect of the physical thickness of the cluster. We have 
assigned all members a distance of 133 pc, although there will naturally be some variation. However, 
this effect is also relatively small. From Section 4.1, the volumetric number density falls off with 
radius approximately as a King model with core and tidal radii of 2.1 and 19 pc, respectively. 
Consider the front half of a spherical cluster with such a density distribution. It may readily be 
shown that the mean distance from the plane of the sky of any cluster member is d = 2.6 pc. For 
a cluster at mean distance D, the induced magnitude spread is 5 log [{D + d)/D], which is 0.04 in 
our case. Although the actual spread in magnitudes is not Gaussian, we have added this figure in 
quadrature to both aj and ax, and rerun the analysis. Again, the binarity is unaffected. 

The errors due to both photometry and finite cluster thickness induce a symmetric spread in 
stellar magnitudes. That is, they scatter as many sources below the fiducial isochrone as above it. 
Thus, they cannot reduce the estimated binarity, which stems from an excess of stars above the 
isochrone. One systematic error that would affect b is an overestimation of the cluster distance. 
If D were lowered, the absolute magnitudes of all sources would decrease equally, as would the 
inferred 5-value. Quantitatively, the distance would have to decrease by about 15 pc to bring the 
binary fraction down to the field-star result for G-dwarf primaries. An error of this size for the 
avera ge distance is excluded by current observations, for which the estimated uncertainty is only 
1 pc (ISoderblom et al.lbond ). 



Since our method relies solely on photometry to assess binarity, we cannot distinguish between 
physically linked pairs and chance alignments. As mentioned in Section 3.4, the resolution limit of 
our data is 10", or Aro = 1400 AU at the distance of the Pleiades. Consider a star at a radius r 
from the cluster center. Its average number of neighbors within Aro is vr Ar^ n<j(r), where ns{r) is 
the projected surface number density of the cluster. Since each ring of width dr contains 2TTnsr dr 
stars, integration over all members yields the total number of chance alignments: 



iVchanco = 2 7r2Ar2 / nl{r)rdr . 
Jo 

where R is the cluster's outer radius. Using ns{r) from Figure 10, we findA'^chance 
chance alignments have no quantitative impact. 



(44) 
2.4. Thus, 



Yet another source of systematic e rror i s the cluster age. We have adopted the lithium-based 
figure of 125 Myr from Istauffer et al.l ( 1998). Earlier est i mate s, using the main-sequence turnoff, 
yielded a range of answers. For example, iMeynet et al.l (119931 ) found 100 Myr. Even this minor 
reduction affects our results, since it lifts the low-mass end of the isochrone toward higher luminosity. 
For an age of 100 Myr, our analysis gives 6 = 0.57 ± 0.02 and c = 0.28 ± 0.06. The binary fract ion 
is augmented to 0.64 when we include resolved pairs. From Figure 3 of IStauffer et al.l (|l998l ). a 
100 Myr age corresponds to a lithium edge at Mj = 11.7, or / = 17.3 at the Pleiades distance. 



Such a result seems incompatible with the lithium data shown in Figure 2 of lStauffer et al.l (|l998l ) 
but the total number of observations is relatively small. 
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We conclude that the enhanced binarity is a real effect. What this fact tells us about the 
origin of open clusters remains to be seen. Our next step in addressing this basic issue is to 
try and account theoretically for the empirical properties just obtained through our statistical 
analysis. We will ascertain, using direct numerical simulations, the range of initial states that 
can relax dynamically to the present-day Pleiades. Such a study will bring us one step closer to 
understanding the molecular cloud environments that give rise to open clusters. 

We are grateful to Eric Huff, who first suggested the use of the maximum likelihood technique 
for this problem and provided continued insight. Wc also benefited from conversations with James 
Graham and Geoff Marcy. The referee's comments helped to improve the final presentation. This 
research was supported by NSF grant AST-0639743. 
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A. Interpreting the Gini Coefficient 

In Section 4.2, we introduced the Gini coefficient geometrically, as twice the area between the 
Im — In curve and the diagnonal line representing zero mass segregation. Alternatively, G may be 
defined in terms of the mean mass difference between radial shells in the cluster. Here we describe 
more precisely, and prove the equivalence of, this second interpretation. 

Altering our previous notation, we now let m{r) be the average system mass (primary plus 
secondary) in a shell with outer radius r. If there are many shells, then we may define a system 
number density n{r), such that the number of systems between r and r + dr is n(r) dr. The total 
number of systems at all radii is 



/•oo 

= / n{r) dr . (Al) 
Jo 



Ntot 

'0 

This total was called /xtot in the text. The average system mass throughout the entire cluster is 



Nt, 



1 f°° 

- — / m(r) n(r) dr . (A2) 
not Jo 

We will be concerned with the relative mean difference in the mass of shells. This is 

-j^ /-ex: poo 

A = „ / dr dr' \m(r) — m(r')\ nir) n(r') . (-^-3) 
iTiNtotJo Jo 

We will prove that the Gini coefficient, as defined in Section 4.2, is also A/2. 

In the geometric definition of G, we utilized the cumulative fractional number /^r and cumu- 
lative fractional mass /m- These may be written in terms of the system number density: 

fN{r) = ^ rn{r')dr' , (A4) 

fM{r) = _ / m{r')n{r')dr' . (A5) 
m iVtot Jo 

We will later need the differentials of /;v and /m, which are 

dfN = "^dr, (A6) 

-'Vtot 

dfM = "^^^ dr . (A7) 
In terms of /at and /m, the geometric definition of G is then 

G = 2 f\fM - In) dfN . (A8) 
Jo 

If the mass is centrally concentrated, as expected in a real stellar cluster, then /m > /n at all 
radii, and G > 0. Hypothetical clusters in which larger masses are preferentially located farther 
from the center would have G < 0. 
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We now return to equation (A3) and manipulate it to obtain G, as defined above. A central 

assumption wc will make is that m(r) declines monotonically with r. We may then expand the 
righthand side of equation (A3). We split the integral over r' into two parts, one for r' < r and the 
other for r' > r. Under our assumption, m{r') > m{r) in the first integral, and m{r') < m{r) in 
the second. After pulling n(r) from the r'-integration, we further split the integrands to find 



A = 



1 r°° 

— / (-Xi + X2 + X3 - X4) n{r) dr . (A9) 
'Vtot Jo 



The first two terms of the integrand are 

m{r) 



Xi = ^ / n{r')dT' = m(r)/iv(r) , (AlO) 







The third term is 



I2 = — — / m{r')n{r') dr' = mfMi'r) ■ (-^H) 

Ntot Jo 



X3 = !^ / n{r')dr' (A12) 



A'tot Jr 

^ rn{r')dr' - ^ rn{r')dr' (A13) 

Jvtot Jo JVtot Jo 

m(r) — m{r) /^{r) , (^14) 



while the fourth is 



1 f°° 

X4 = -— / m{r')n{r')dr' (A15) 

A'tot Jr 

= — — / m(r') ra(r') dr' — — - — / m(r') n{r') dr' (-^16) 
J^tot Jo -'not Jo 

= fh - fhfMir) . (A17) 

Putting equations (AlO), (All), (A14), and (A17) back into equation (A9) yields 

A = _ — / [2m/M(r) — 2m(r)/Ar(r) + m(r) — m] n{r)dr (^18) 
mNtot Jo 

= 2 (^J^ fMdfN - ^ fNdfM^ , (A19) 

where we have used both equations (Al) and (A2) to cancel the last two terms on the right side of 
equation (A18), and equations (A6) and (A7) to transform the remaining two terms. The second 
integral in the last equation is 

/ fNdfM = / difNfM) - / fMdfN (A20) 

Jo Jfj^fM=0 Jo 



JO 



fMdfN, (A21) 
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since /m and /at attain their upper and lower bounds simultaneously. Using this result, equa- 
tion (A19) becomes 

A = 2 f 2 / fMdfN - l] . (A22) 

Finally, we note that 







Thus, equation (A22) becomes 



as claimed. 



Jo 



2 / fNdfN = l. (A23) 



A = 4 (^l^'fMdfN - J^'fNdfN^ (A24) 

= 4 / (/m - /at) dfN (A25) 
Jo 

= 2G , (A26) 
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Fig. 1. — Near-infrared color-magnitude diagram for the Pleiades. Small dots represent the 1245 

stars in our sample. Open circles are the 41 likely sub-stellar objects which have been removed from 
the sample. Filled circles arc the 11 brightest stars, which arc likely post-main-sequence objects. 
The 125 Myr isochrone for stars with masses between 0.08 and 4.0 Mq is shown as the smooth, 
solid curve. 
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Fig. 2. — Observational error in t he i^T-band measurem ents as a function of absolute magnitude 
for all 1417 stars in the catalog of lStauffer et al.l (j2007l ). The smooth curve is the approximate fit 
used in our maximum likelihood analysis. 
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Stellar Mass log m 

Fig. 3. — Sample synthetic data results for the single star mass function (j)[ra). In both panels, the 
smooth curve is the input function from which the synthetic data were drawn. Shown also are the 
best-fit values, along with errors, for our discrete representation of the function. In the left panel, 
the input ^(m) is a power law with slope -2.8; in the right panel, it is a log-normal function with 
peak mo = 0.2 and width = 0.4. 
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Fig. 4. — Weighted sum of the biases as a function of the regularization parameter A. The synthetic 
input here was a log-normal function with 1245 sources. The dashed, horizontal line is drawn at 
xl = N, where A/" = 21 is the number of free parameters in our fitting. The short, vertical arrow 
indicates the final value of A used for this synthetic dataset. 
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Fig. 5. — Comparison of our fitted correlation coefficient c with 7, the imposed degree of correlation 
in the synthetic dataset. The gray area indicates the region in which the Pleiades most likely falls. 
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Fig. 6. — Systematic errors in parameters of the log-normal fit to (i){m), plotted as a function of 
the synthetic binary correlation 7. Filled circles show Anio, the error in ruo- Open circles show 
Aajn, the error in am- The gray areas indicate the regions in which the Pleiades most likely falls. 
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Fig. 7. — Best-fit single star probability density ^(m) for the Pleiades. Actual bin values yj/Amj 
are shown with associated errors. The smooth curve is a log-normal approximation to the results. 
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Fig. 8. — Comparison of log-normal fits to the primary probability density log(j)p{mp) (solid curve) 
and the single star probability density log (j){m) (dashed curve). The primary function peaks at 
larger mass and has a smaller width. Note that 4>p{m) includes single stars. 
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Fig. 9. — Comparison of the Pleiades single star probability density to the field-star initial mass 
function of iKroupal (j200lh . where the latter has been shifted upward for clarity . Shown here is the 
log-normal approximation to log (j){m) (solid curve), augmented with the data of lBihain et all (|2006l ) 
for low-mass members and brown dwarfs (open circles) and our 11 brightest, post-main-sequence 
stars (filled circles). 
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Fig. 10. — Surface density distribution in the Pleiades. The filled circles represent the surface 
number density (pc^^), displayed on a logarithmic scale. Open circles are the mass density, in 
Mq pc~^. The solid and dashed smooth curves are King model fits. 
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Fig. 11. — Volume density profiles. The filled circles represent the number density (pc~^), again 
displayed logarithmically. Open circles are the mass density, in Mq pc~^. The solid and dashed 
smooth curves are the same King model fits as in Figure 10, but now deprojected into three- 
dimensional space. 



-39- 




4 6 8 10 
Projected Radius (pc) 



12 



Fig. 12. — Average system mass (primaries plus secondaries) as a function of projected cluster 
radius. The bins here have constant radial width. 



-40- 




Fig. 13. — Fractional mass versus fractional number for the Pleiades. The data poi n ts an d error 
bars, along with the solid curve, utilize all sources in the catalog of IStauffer et al.l (j2007l ). The 
dashed curve shows the result when the 11 brightest stars are removed. In both cases, the radial 
bins contain roughly equal numbers of stars. Finally, the dotted diagonal is the hypothetical result 
for no mass segregation. 
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Fig. 14. — Fractional mass versus fract ional number for the ONC. The sohd curve was computed 
using all the ONC sources analyzed bv lHuff &: Stahleii (j2006l ). The dashed curve shows the result 
if only the 4 Trapezium stars are removed. As in Figure 13, the radial bins contain roughly equal 
numbers of stars. The dotted diagonal again shows the hypthetical condition of no mass segregation. 



